clear all
set more off

global data 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_dta"
global figures 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_log"
global in 		"R:\SharedProjects\Shared2020-070\2016\input"


cap log close
log using $figures\M_file,replace t

**************************************************************************************************************
cap program drop pdm1m0
program pdm1m0
	args x1 x2 r
	g pDi=binormal(`x1',`x2',`r')/normal(`x2')						/*This is Pr(D=1|A=1)*/
	g m1i=normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',`x1',`r') + ///
		`r'*normalden(`x2')*(normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',`x1',`r')
	gen m0i=-normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',-`x1',-`r') + ///
		`r'*normalden(`x2')*(1-normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',-`x1',-`r')	
	scalar sigma=sqrt(1+s_csi^2^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)

	qui su part1 if F==1 & A==1
	scalar p1_1=r(mean)
	qui su part1 if F==0 & A==1
	scalar p1_0=r(mean)
	scalar part1m=p1_1-p1_0
	qui su part2 if F==1 & A==1
	scalar p2_1=r(mean)
	qui su part2 if F==0 & A==1
	scalar p2_0=r(mean)
	scalar part2m=p2_1-p2_0
	qui su part3 if F==1 & A==1
	scalar part3m=r(mean)
	scalar drop sigma p1_1 p1_0 p2_1 p2_0
	drop pDi m1i m0i part1 part2 part3
end
**************************************************************************************************************

*********************************************************************************************************************************
/*CORRECTION FACTORS*/
cap program drop get_correction
program get_correction
	scalar cfl   =bF_L-((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2))) 
	scalar cfa   =bF_A-((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))
	scalar cfr   =bF_R-delta
	scalar cfcla =rho_LA-(sqrt(1+s_omega^2+s_psi^2)/sqrt(1+s_omega^2+s_psi^2+s_v^2))
	scalar cfclr =rho_LR-(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2)))
	scalar cfcra =rho_RA-(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2+s_v^2)))
end
***********************************************************************************************************************************

cap program drop get_baseline_pars
program get_baseline_pars
	preserve
		u $data\R_estimates_DWMD,clear
		destring DWMD,gen(beta)
		replace namev="alpha" if namev=="alpha1"
		ren namev namepar
		*********************************************************************************************************************************
		*BASELINE PARAMETERS
		*********************************************************************************************************************************
		su beta if namepar=="pi"
		scalar mu		=r(mean) 
		su beta if namepar=="gamma"
		scalar gamma	=r(mean)
		su beta if namepar=="tau"
		scalar tau		=r(mean)
		su beta if namepar=="sigma_omega"
		scalar s_omega	=r(mean)
		su beta if namepar=="sigma_psi"
		scalar s_psi	=r(mean)
		su beta if namepar=="sigma_v"
		scalar s_v		=r(mean)
		su beta if namepar=="sigma_csi"
		scalar s_csi	=r(mean)
		su beta if namepar=="alpha"
		scalar alpha1	=r(mean)
		scalar sigma	=sqrt(1+s_csi^2)			/*This is [1+var(noise_SSA)] */	
		***********************************************************************************************************************************
		*********************************************************************************************************************************
	restore
end


cd $data
***********************************************************************************************************************************************************************************
u $data\after_ML, clear

ren female F
gen Aw=1-R
ren logbs Y

ren a1 Lhat										/*This includes the effect of gender*/
ren a2 Ahat										/*This includes the effect of gender*/		
ren a3 Rhat										/*This includes the effect of gender*/

ren a1_institutional_year 		xbL_inst
ren a2_institutional_year 		xbA_inst
ren a3_institutional_year 		xbR_inst
gen xbL_noninst 	=Lhat-xbL_inst-bF_L*F		/*This excludes the effect of gender to avoid double counting in the exercise below*/
gen xbA_noninst 	=Ahat-xbA_inst-bF_A*F		/*This excludes the effect of gender to avoid double counting in the exercise below*/
gen xbR_noninst 	=Rhat-xbR_inst-bF_R*F		/*This excludes the effect of gender to avoid double counting in the exercise below*/

ren a1_demo 	xbL_demo						/*This includes the effect of gender*/
replace xbL_demo=xbL_demo-bF_L*F				/*This excludes it to avoid double counting*/
ren a2_demo 	xbA_demo						/*This includes the effect of gender*/
replace xbA_demo=xbA_demo-bF_A*F				/*This excludes it to avoid double counting*/
ren a3_demo 	xbR_demo						/*This includes the effect of gender*/
replace xbR_demo=xbR_demo-bF_R*F				/*This excludes it to avoid double counting*/

ren a1_labor 	xbL_labor
ren a1_health 	xbL_health
ren a2_labor 	xbA_labor
ren a2_health 	xbA_health
ren a3_labor 	xbR_labor
ren a3_health 	xbR_health

scalar rho_LA=rho_A_L
scalar rho_RA=rho_A_R
scalar rho_LR=rho_L_R
scalar betaL=bF_L
scalar betaA=bF_A
scalar betaR=bF_R

g Khat=-Lhat

**********************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global spec_L 	F college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm 
global spec_R 	F college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx*
global spec_A 	F college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm Y 
***********************************************************************************************************************************************************************************

get_baseline_pars

gen xbL=Lhat-bF_L*F			/*subtract the effect of gender because it is going to add it "structurally" later*/
gen xbA=Ahat-bF_A*F			/*subtract the effect of gender because it is going to add it "structurally" later*/
gen xbR=Rhat-bF_R*F			/*subtract the effect of gender because it is going to add it "structurally" later*/	

gen Dhat=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

pdm1m0 Dhat Ahat rho_DA

scalar delta=part1m+part2m+part3m*alpha1 	

get_correction

gen cfid=1
replace cfid=sum(cfid)
				
matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
matrix W2=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat Ahat Rhat) if !missing(Rhat), replace
putmata cfid U2=(Khat Ahat Rhat) if !missing(Rhat), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
	st_matrix("Z",Z)
	st_matrix("K",K)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

*********************************************************************************************************************************
*****Recall: Pr(R=1|A=1)=Pr(R=1|A=1,L=1)Pr(L=1|A=1)+Pr(R=1|A=1,L=0)Pr(L=0|A=1)
g P_R_given_A2=binormal(Rhat,Ahat,rho_RA)/normal(Ahat)			/*This is Pr(R=1|A=1) */
g P_L_given_A2=binormal(Lhat,Ahat,rho_LA)/normal(Ahat)			/*This is Pr(L=1|A=1) */
g P_notL_given_A2=binormal(-Lhat,Ahat,-rho_LA)/normal(Ahat)		/*This is Pr(L=0|A=1) */
g P_R_given_AL2=Z/binormal(Ahat,Lhat,rho_LA)					/*This is Pr(R=1|L=1,A=1)*/
g P_R_given_AnotL2=K/binormal(Ahat,-Lhat,-rho_LA)				/*This is Pr(R=0|L=0,A=1)*/
*********************************************************************************************************************************

*********************************************************************************************************************************
								/*CHECK THAT IT WORKS*/
g lhs=P_R_given_A2
g rhs=P_R_given_AL2*P_L_given_A2+P_R_given_AnotL2*P_notL_given_A2
su lhs rhs						/*it does*/
*********************************************************************************************************************************

*********************************************************************************************************************************
g Type1=P_R_given_AL			/*This is Pr(type I error)*/
g Type2=1-P_R_given_AnotL		/*This is Pr(type II error)*/

su Type1 if F==1 & A==1			
su Type1 if F==0 & A==1
su Type2 if F==1 & A==1
su Type2 if F==0 & A==1
*********************************************************************************************************************************

*********************************************************************************************************************************
drop Z K
*********************************************************************************************************************************

save $data\cfactual,replace


u $data\cfactual,clear	
*********************************************************************************************************************************
****C O U N T E R F A C T U A L S
*********************************************************************************************************************************
*********************************************************************************************************************************

*(BAN THE BOX EXPERIMENT)
*********************************************************************************************************************************
*********************************************************************************************************************************
*BASELINE PARAMETERS
*********************************************************************************************************************************
get_baseline_pars
***********************************************************************************************************************************
***********************************************************************************************************************************		

probit F college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx*
predict F_P,pr

g Lhat_Wsex=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
g Ahat_Wsex=xbA+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
g Dhat_Wsex=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F_P
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

pdm1m0 Dhat_Wsex Ahat_Wsex rho_DA
scalar delta=part1m+part2m+part3m*alpha1 		

g Rhat_Wsex=xbR+(delta+cfr)*F_P												/*delta already includes the scaling by sigma*/
g Khat_Wsex=-Lhat_Wsex

matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)			/*The correlation coefficients don't change, so keep baseline*/
matrix W2=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)

putmata cfid U1=(Lhat_Wsex Ahat_Wsex Rhat_Wsex) if !missing(Rhat_Wsex), replace
putmata cfid U2=(Khat_Wsex Ahat_Wsex Rhat_Wsex) if !missing(Rhat_Wsex), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

g Type1_Wsex=Z   /binormal(Ahat_Wsex,Lhat_Wsex,rho_LA)
g Type2_Wsex=1-(K/binormal(Ahat_Wsex,Khat_Wsex,-rho_LA))

su Type1* if F==1 & A==1
su Type1* if F==0 & A==1
su Type1* if A==1



*********************************************************************************************************************************
* PI=GAMMA=TAU=0 (only taste discrimination)
*********************************************************************************************************************************
*********************************************************************************************************************************
*BASELINE PARAMETERS
*********************************************************************************************************************************
get_baseline_pars	
***********************************************************************************************************************************
***********************************************************************************************************************************		
scalar tau=0
scalar mu=0
scalar gamma=0


g Lhat_pgt=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
g Ahat_pgt=xbA+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
g Dhat_pgt=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

pdm1m0 Dhat_pgt Ahat_pgt rho_DA
scalar delta=part1m+part2m+part3m*alpha1 		

g Rhat_pgt=xbR+(delta+cfr)*F
g Khat_pgt=-Lhat_pgt


matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
matrix W2=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat_pgt Ahat_pgt Rhat_pgt) if !missing(Rhat_pgt), replace
putmata cfid U2=(Khat_pgt Ahat_pgt Rhat_pgt) if !missing(Rhat_pgt), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

	
g Type1_pgt=Z   /binormal(Ahat_pgt,Lhat_pgt,rho_LA)
g Type2_pgt=1-(K/binormal(Ahat_pgt,Khat_pgt,-rho_LA))

drop Z	K Lhat_* Ahat_* Rhat_* Khat_*

su Type* if F==1 & A==1
su Type* if F==0 & A==1


*********************************************************************************************************************************
* PI=GAMMA=TAU=ALPHA1=0 (only X-heterogeneity)
*********************************************************************************************************************************
*********************************************************************************************************************************
*BASELINE PARAMETERS
*********************************************************************************************************************************
get_baseline_pars	
***********************************************************************************************************************************
***********************************************************************************************************************************		
scalar tau=0
scalar mu=0
scalar gamma=0
scalar alpha1=0


g Lhat_pgta=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
g Ahat_pgta=xbA+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
g Dhat_pgta=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

pdm1m0 Dhat_pgta Ahat_pgta rho_DA
scalar delta=part1m+part2m+part3m*alpha1 		

g Rhat_pgta=xbR+(delta+cfr)*F
g Khat_pgta=-Lhat_pgta


matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
matrix W2=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat_pgta Ahat_pgta Rhat_pgta) if !missing(Rhat_pgta), replace
putmata cfid U2=(Khat_pgta Ahat_pgta Rhat_pgta) if !missing(Rhat_pgta), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

	
g Type1_pgta=Z   /binormal(Ahat_pgta,Lhat_pgta,rho_LA)
g Type2_pgta=1-(K/binormal(Ahat_pgta,Khat_pgta,-rho_LA))

drop Z	K Lhat_* Ahat_* Rhat_* Khat_*

su Type* if F==1 & A==1
su Type* if F==0 & A==1



*********************************************************************************************************************************
* ALPHA1=0 (only statistical discrimination)
*********************************************************************************************************************************
*********************************************************************************************************************************
*BASELINE PARAMETERS
*********************************************************************************************************************************
get_baseline_pars	
***********************************************************************************************************************************
***********************************************************************************************************************************		
scalar alpha1=0

g Lhat_alpha=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
g Ahat_alpha=xbA+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
g Dhat_alpha=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

pdm1m0 Dhat_alpha Ahat_alpha rho_DA
scalar delta=part1m+part2m+part3m*alpha1 		
				
g Rhat_alpha=xbR+(delta+cfr)*F
g Khat_alpha=-Lhat_alpha

matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
matrix W2=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat_alpha Ahat_alpha Rhat_alpha) if !missing(Rhat_alpha), replace
putmata cfid U2=(Khat_alpha Ahat_alpha Rhat_alpha) if !missing(Rhat_alpha), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

	
g Type1_alpha=Z   /binormal(Ahat_alpha,Lhat_alpha,rho_LA)
g Type2_alpha=1-(K/binormal(Ahat_alpha,Khat_alpha,-rho_LA))

drop Z	K Lhat_* Ahat_* Rhat_* Khat_*

su Type1* if F==1 & A==1
su Type1* if F==0 & A==1


******
******
******
******

get_baseline_pars	
scalar alpha1=0
scalar mu=0
scalar gamma=0
scalar tau=0

su xbL_l
g xlml=r(mean)
su xbL_d
g xlmd=r(mean)
su xbL_h
g xlmh=r(mean)
su xbA_l
g xaml=r(mean)
su xbA_d
g xamd=r(mean)
su xbA_h
g xamh=r(mean)
su xbR_l
g xrml=r(mean)
su xbR_d
g xrmd=r(mean)
su xbR_h
g xrmh=r(mean)

g Lhat_all0_exp1=xlml+xlmd+xbL_health+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2))))*F
replace Lhat_all0_exp1=. if Lhat==.
g Ahat_all0_exp1=xaml+xamd+xbA_health+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2))))*F
replace Ahat_all0_exp1=. if Ahat==.
g Dhat_all0_exp1=Lhat_all0_exp1*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
pdm1m0 Dhat_all0_exp1 Ahat_all0_exp1 rho_DA
scalar delta=part1m+part2m+part3m*alpha1 		
g Rhat_all0_exp1=xrml+xrmd+xbR_health+delta*F
replace Rhat_all0_exp1=. if Rhat==.

matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat_all0_exp1 Ahat_all0_exp1 Rhat_all0_exp1) if !missing(Rhat_all0_exp1), replace
mata
	R1=vech(st_matrix("W1"))'
	Z=mvnormal(U1,R1)
end
getmata Z, id(cfid) replace
g Type1_all0_exp1=Z   /binormal(Ahat_all0_exp1,Lhat_all0_exp1,rho_LA)
drop Z	Lhat_* Ahat_* Rhat_*  
su Type1* if F==1 & A==1
su Type1* if F==0 & A==1

g Lhat_all0_exp2=xlmd+xbL_labor+xbL_health+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2))))*F
replace Lhat_all0_exp2=. if Lhat==.
g Ahat_all0_exp2=xamd+xbA_labor+xbA_health+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2))))*F
replace Ahat_all0_exp2=. if Ahat==.
g Dhat_all0_exp2=Lhat_all0_exp2*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
pdm1m0 Dhat_all0_exp2 Ahat_all0_exp2 rho_DA
scalar delta=part1m+part2m+part3m*alpha1 		
g Rhat_all0_exp2=xrmd+xbR_labor+xbR_health+delta*F
replace Rhat_all0_exp2=. if Rhat==.

matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat_all0_exp2 Ahat_all0_exp2 Rhat_all0_exp2) if !missing(Rhat_all0_exp2), replace
mata
	R1=vech(st_matrix("W1"))'
	Z=mvnormal(U1,R1)
end
getmata Z, id(cfid) replace
g Type1_all0_exp2=Z   /binormal(Ahat_all0_exp2,Lhat_all0_exp2,rho_LA)
drop Z	Lhat_* Ahat_* Rhat_*  
su Type1* if F==1 & A==1
su Type1* if F==0 & A==1

g xbL_demo_noage=xbL_demo-a1_age
g xbA_demo_noage=xbA_demo-a2_age
g xbR_demo_noage=xbR_demo-a3_age
su xbL_demo_noage
replace xlmd=r(mean)
su xbA_demo_noage
replace xamd=r(mean)
su xbR_demo_noage
replace xrmd=r(mean)

g Lhat_all0_exp3=xlmd+a1_age+xbL_labor+xbL_health+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2))))*F
replace Lhat_all0_exp3=. if Lhat==.
g Ahat_all0_exp3=xamd+a2_age+xbA_labor+xbA_health+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2))))*F
replace Ahat_all0_exp3=. if Ahat==.
g Dhat_all0_exp3=Lhat_all0_exp3*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
pdm1m0 Dhat_all0_exp3 Ahat_all0_exp3 rho_DA
scalar delta=part1m+part2m+part3m*alpha1 		
g Rhat_all0_exp3=xrmd+a3_age+xbR_labor+xbR_health+delta*F
replace Rhat_all0_exp3=. if Rhat==.

matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat_all0_exp3 Ahat_all0_exp3 Rhat_all0_exp3) if !missing(Rhat_all0_exp3), replace
mata
	R1=vech(st_matrix("W1"))'
	Z=mvnormal(U1,R1)
end
getmata Z, id(cfid) replace
g Type1_all0_exp3=Z   /binormal(Ahat_all0_exp3,Lhat_all0_exp3,rho_LA)	
drop Z	Lhat_* Ahat_* Rhat_*  
su Type1* if F==1 & A==1				/*This is the same as "institutional variables"*/
su Type1* if F==0 & A==1


*********************************************************************************************************************************
* Lbar_i=Lbar_SSA
*********************************************************************************************************************************
*********************************************************************************************************************************
*BASELINE PARAMETERS
*********************************************************************************************************************************
get_baseline_pars	
***********************************************************************************************************************************
***********************************************************************************************************************************		
scalar s_psi_new	=0
scalar gamma		=0

scalar rho_LA_new	=(sqrt(1+s_omega^2+s_psi_new^2)/sqrt(1+s_omega^2+s_psi_new^2+s_v^2))+cfcla
scalar rho_LR_new	=(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi_new^2)))+cfclr
scalar rho_RA_new	=(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi_new^2+s_v^2)))+cfcra

scalar sL		=sqrt(1+s_omega^2+s_psi^2)
scalar sL_new	=sqrt(1+s_omega^2+s_psi_new^2)
scalar sA		=sqrt(1+s_omega^2+s_psi^2+s_v^2)
scalar sA_new	=sqrt(1+s_omega^2+s_psi_new^2+s_v^2)
g Lhat_ssa=Lhat*(sL/sL_new)
g Ahat_ssa=Ahat*(sA/sA_new)
g Dhat_ssa=xbL*sL+mu*F
scalar rho_DA=1/sA_new

pdm1m0 Dhat_ssa Ahat_ssa rho_DA
		
scalar delta=part1m+part2m+part3m*alpha1 		

g Rhat_ssa=xbR+(delta+cfr)*F
g Khat_ssa=-Lhat_ssa

matrix W1=(1,  rho_LA_new,  rho_LR_new \  rho_LA_new,1,rho_RA_new \  rho_LR_new, rho_RA_new, 1)
matrix W2=(1, -rho_LA_new, -rho_LR_new \ -rho_LA_new,1,rho_RA_new \ -rho_LR_new, rho_RA_new, 1)
putmata cfid U1=(Lhat_ssa Ahat_ssa Rhat_ssa) if !missing(Rhat_ssa), replace
putmata cfid U2=(Khat_ssa Ahat_ssa Rhat_ssa) if !missing(Rhat_ssa), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

	
g Type1_ssa=Z   /binormal(Ahat_ssa,Lhat_ssa,rho_LA_new)
g Type2_ssa=1-(K/binormal(Ahat_ssa,Khat_ssa,-rho_LA_new))

drop Z	K Lhat_* Ahat_* Rhat_* Khat_*


*********************************************************************************************************************************
* S_OMEGA=0
*********************************************************************************************************************************
*********************************************************************************************************************************
*BASELINE PARAMETERS
*********************************************************************************************************************************
get_baseline_pars	
***********************************************************************************************************************************
***********************************************************************************************************************************		
scalar s_omega_new	=0

scalar rho_LA_new	=(sqrt(1+s_omega_new^2+s_psi^2)/sqrt(1+s_omega_new^2+s_psi^2+s_v^2))+cfcla
scalar rho_LR_new	=(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega_new^2+s_psi^2)))+cfclr
scalar rho_RA_new	=(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega_new^2+s_psi^2+s_v^2)))+cfcra

scalar sL		=sqrt(1+s_omega^2+s_psi^2)
scalar sL_new	=sqrt(1+s_omega_new^2+s_psi^2)
scalar sA		=sqrt(1+s_omega^2+s_psi^2+s_v^2)
scalar sA_new	=sqrt(1+s_omega_new^2+s_psi^2+s_v^2)
g Lhat_omega=Lhat*(sL/sL_new)
g Ahat_omega=Ahat*(sA/sA_new)

g Dhat_omega=xbL*sL+mu*F
scalar rho_DA=1/sA_new

pdm1m0 Dhat_omega Ahat_omega rho_DA
		
scalar delta=part1m+part2m+part3m*alpha1 		

g Rhat_omega=xbR+(delta+cfr)*F
g Khat_omega=-Lhat_omega

matrix W1=(1,  rho_LA_new,  rho_LR_new \  rho_LA_new,1,rho_RA_new \  rho_LR_new, rho_RA_new, 1)
matrix W2=(1, -rho_LA_new, -rho_LR_new \ -rho_LA_new,1,rho_RA_new \ -rho_LR_new, rho_RA_new, 1)
putmata cfid U1=(Lhat_omega Ahat_omega Rhat_omega) if !missing(Rhat_omega), replace
putmata cfid U2=(Khat_omega Ahat_omega Rhat_omega) if !missing(Rhat_omega), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

	
g Type1_omega=Z   /binormal(Ahat_omega,Lhat_omega,rho_LA_new)
g Type2_omega=1-(K/binormal(Ahat_omega,Khat_omega,-rho_LA_new))

drop Z	K Lhat_* Ahat_* Rhat_* Khat_*


*********************************************************************************************************************************
* Lbar_i=Lbar_SSA, gamma=0
*********************************************************************************************************************************
*********************************************************************************************************************************
*BASELINE PARAMETERS
*********************************************************************************************************************************
get_baseline_pars	
***********************************************************************************************************************************
***********************************************************************************************************************************		
scalar s_psi_new	=0
scalar gamma		=0
scalar s_omega_new	=0


scalar rho_LA_new	=(sqrt(1+s_omega_new^2+s_psi_new^2)/sqrt(1+s_omega_new^2+s_psi_new^2+s_v^2))+cfcla
scalar rho_LR_new	=(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega_new^2+s_psi_new^2)))+cfclr
scalar rho_RA_new	=(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega_new^2+s_psi_new^2+s_v^2)))+cfcra

scalar sL		=sqrt(1+s_omega^2+s_psi^2)
scalar sL_new	=sqrt(1+s_omega_new^2+s_psi_new^2)
scalar sA		=sqrt(1+s_omega^2+s_psi^2+s_v^2)
scalar sA_new	=sqrt(1+s_omega_new^2+s_psi_new^2+s_v^2)
g Lhat_ssa2=Lhat*(sL/sL_new)
g Ahat_ssa2=Ahat*(sA/sA_new)

g Dhat_ssa2=xbL*sL+mu*F
scalar rho_DA=1/sA_new

pdm1m0 Dhat_ssa2 Ahat_ssa2 rho_DA
		
scalar delta=part1m+part2m+part3m*alpha1 		

g Rhat_ssa2=xbR+(delta+cfr)*F			/*The standard deviation of R doesn't depend on s_omega*/
g Khat_ssa2=-Lhat_ssa2

matrix W1=(1,  rho_LA_new,  rho_LR_new \  rho_LA_new,1,rho_RA_new \  rho_LR_new, rho_RA_new, 1)
matrix W2=(1, -rho_LA_new, -rho_LR_new \ -rho_LA_new,1,rho_RA_new \ -rho_LR_new, rho_RA_new, 1)
putmata cfid U1=(Lhat_ssa2 Ahat_ssa2 Rhat_ssa2) if !missing(Rhat_ssa2), replace
putmata cfid U2=(Khat_ssa2 Ahat_ssa2 Rhat_ssa2) if !missing(Rhat_ssa2), replace

mata
	R1=vech(st_matrix("W1"))'
	R2=vech(st_matrix("W2"))'
	Z=mvnormal(U1,R1)
	K=mvnormal(U2,R2)
end

getmata Z, id(cfid) replace
getmata K, id(cfid) replace

	
g Type1_ssa2=Z   /binormal(Ahat_ssa2,Lhat_ssa2,rho_LA_new)
g Type2_ssa2=1-(K/binormal(Ahat_ssa2,Khat_ssa2,-rho_LA_new))

drop Z	K Lhat_* Ahat_* Rhat_* Khat_*

su Type1* if F==1 & A==1
su Type1* if F==0 & A==1

****************************************************************************************************************

preserve
	keep Type* F A
	save cfactual_results,replace
restore


u $data\cfactual_results,clear

keep if Type1!=.

kdensity Type1 if F==1,gen(x1_baseline d1_baseline) nogr
kdensity Type1 if F==0,at(x1_baseline) gen(x0_baseline d0_baseline) nogr

lab var d1_baseline "Women" 
lab var d0_baseline "Men"
lab var x1_baseline "Type I error, baseline"
scatter d1_baseline d0_baseline x1_baseline,c(l l) s(i i) clp(solid dash) legend(pos(6) col(2)) saving(a1,replace) graphr(c(white))

kdensity Type1_all0_exp3 if F==1,gen(x1_inst d1_inst) nogr
kdensity Type1_all0_exp3 if F==0,at(x1_inst) gen(x0_inst d0_inst) nogr
lab var d1_inst "Women" 
lab var d0_inst "Men"
lab var x1_inst "Type I error, only inst.var. heterog."
scatter d1_inst d0_inst x1_inst,c(l l) s(i i) clp(solid dash) legend(pos(6) col(2)) saving(a3,replace) graphr(c(white))

graph combine a3.gph a1.gph,graphr(c(white))
graph export $figures\figure2.pdf,replace




****Create Latex Tables

u $data\cfactual_results,clear
keep if A==1
gen M=1-F
gcollapse (mean) Type1*
gen M=3
save $data\agg_means,replace

u $data\cfactual_results,clear
keep if A==1
gen M=1-F
gcollapse (mean) Type1*,by(M)
set obs 3
replace M=999 in 3

append using $data\agg_means


foreach x in Type1_all0_exp1 Type1_all0_exp2 Type1_all0_exp3 Type1_pgta Type1_alpha Type1_pgt Type1 {
	replace `x'=(`x'[_n-2]-`x'[_n-1]) if M==999
}
lab define M 1 "Male" 0 "Female" 3 "All" 999 "Difference",replace
lab values M M

label var Type1 			"1. Estimated Type I Difference"
label var Type1_all0_exp3 	"2. No Discrimination, only institutional var. heterog."
label var Type1_pgta 		"3. No Discrimination, only X-Heterogeneity"
label var Type1_alpha 		"4. Allowing only for Statistical Discrimination"
label var Type1_pgt 		"5. Allowing only for Taste Discrimination"

eststo s_all: estpost tabstat Type1 Type1_all0_exp3 Type1_pgta Type1_alpha Type1_pgt ///
	 if M==3, statistics(mean) columns(statistics)
eststo s0: estpost tabstat Type1 Type1_all0_exp3 Type1_pgta Type1_alpha Type1_pgt  ///
	 if M==0, statistics(mean) columns(statistics)
eststo s1: estpost tabstat Type1 Type1_all0_exp3 Type1_pgta Type1_alpha Type1_pgt  ///
	 if M==1, statistics(mean) columns(statistics)
eststo sd: estpost tabstat Type1 Type1_all0_exp3 Type1_pgta Type1_alpha Type1_pgt  ///
	 if M==999, statistics(mean) columns(statistics)

esttab s_all s0 s1 sd using $figures\table8.tex,   ///
	 cells("mean(fmt(2))")  nodepvar label unstack mtitle("All" "Women" "Men" "Difference") nonumber noobs ///
	 title("Counterfactuals" \label{tab8}) replace ///
	 postfoot(\hline\hline \end{tabular} { \\ \begin{flushleft} \begin{scriptsize} \textit{Note:} \\ ///
	 Row 1: Baseline. \\ ///
	 Row 2: $\alpha=\gamma=\pi=\tau=0$, only institutional heterogeneity. \\ ///
	 Row 3: $\alpha=\gamma=\pi=\tau=0$, all heterogeneity. \\ ///
	 Row 4: $\alpha=0$. \\ ///
	 Row 5: $\gamma=\pi=\tau=0$. \\ ///
	 \end{scriptsize} \end{flushleft} } \end{table} )
eststo clear





*****************************************************Sex experiments
u $data\cfactual_results,clear
keep if A==1
gen M=1-F
gcollapse (mean) Type1*
gen M=3
save $data\agg_means,replace

u $data\cfactual_results,clear
keep if A==1
gen M=1-F
gcollapse (mean) Type1*,by(M)
set obs 3
replace M=999 in 3

append using $data\agg_means


foreach x in Type1 Type1_alpha Type1_Wsex  {
	replace `x'=(`x'[_n-2]-`x'[_n-1]) if M==999
}
lab define M 1 "Male" 0 "Female" 3 "All" 999 "Difference",replace
lab values M M

label var Type1 			"1. Estimated Type I Difference"
label var Type1_alpha 		"2. No Taste-Based Discrimination"
label var Type1_Wsex 		"3. Ban the (Gender) Box"

eststo clear
eststo s_all: estpost tabstat Type1 Type1_alpha Type1_Wsex ///
	 if M==3, statistics(mean) columns(statistics)
eststo s0: estpost tabstat Type1 Type1_alpha Type1_Wsex  ///
	 if M==0, statistics(mean) columns(statistics)
eststo s1: estpost tabstat Type1 Type1_alpha Type1_Wsex  ///
	 if M==1, statistics(mean) columns(statistics)
eststo sd: estpost tabstat Type1 Type1_alpha Type1_Wsex  ///
	 if M==999, statistics(mean) columns(statistics)
	 
esttab s_all s0 s1 sd using $figures\table11.tex,   ///
	 cells("mean(fmt(2))")  nodepvar label unstack mtitle("All" "Women" "Men" "Difference") nonumber noobs ///
	 title("Counterfactuals" \label{tab8sex}) replace ///
	 postfoot(\hline\hline \end{tabular} { \\ \begin{flushleft} \begin{scriptsize} \textit{Note:} \\ ///
	 Row 1: Baseline. \\ ///
	 Row 2: No Taste-Based Discrimination. \\ ///
	 Row 3: Ban the (Gender) Box. \\ ///
	 \end{scriptsize} \end{flushleft} } \end{table} )
eststo clear


*****************************************************
u $data\cfactual_results,clear
keep if A==1
gen M=1-F
gcollapse (mean) Type1*
gen M=3
save $data\agg_means,replace

u $data\cfactual_results,clear
keep if A==1
gen M=1-F
gcollapse (mean) Type1*,by(M)
set obs 3
replace M=999 in 3

append using $data\agg_means

ren Type1_ssa2 Type1_ssa_omega
foreach x in Type1 Type1_ssa Type1_omega Type1_ssa_omega {
	replace `x'=(`x'[_n-2]-`x'[_n-1]) if M==999
}
lab define M 1 "Male" 0 "Female" 3 "All" 999 "Difference",replace
lab values M M

label var Type1 "Estimated Type I Difference"
label var Type1_ssa  "No Heterogeneity in Self-Report Thresholds"
label var Type1_omega "No Noise in Self-Reports"
label var Type1_ssa_omega "No Heterogeneity in Self-Report Thresholds and No Noise in Self-Reports"

eststo clear
eststo s_all: estpost tabstat Type1 Type1_ssa Type1_omega  Type1_ssa_omega ///
	 if M==3, statistics(mean) columns(statistics)
eststo s0: estpost tabstat Type1 Type1_ssa Type1_omega  Type1_ssa_omega  ///
	 if M==0, statistics(mean) columns(statistics)
eststo s1: estpost tabstat Type1 Type1_ssa Type1_omega  Type1_ssa_omega  ///
	 if M==1, statistics(mean) columns(statistics)
eststo sd: estpost tabstat Type1 Type1_ssa Type1_omega  Type1_ssa_omega  ///
	 if M==999, statistics(mean) columns(statistics)

esttab s_all s0 s1 sd using $figures\table9.tex,   ///
	 cells("mean(fmt(2))")  nodepvar label unstack mtitle("All" "Women" "Men" "Difference") nonumber noobs ///
	 title("Counterfactuals" \label{tab8sex}) replace ///
	 postfoot(\hline\hline \end{tabular} { \\ \begin{flushleft} \begin{scriptsize} \textit{Note:} \\ ///
	 Row 1: Baseline. \\ ///
	 Row 2: No Heterogeneity in Self-Report Thresholds. \\ ///
	 Row 3: No Noise in Self-Reports. \\ ///
	 Row 4: No Heterogeneity in Self-Report Thresholds and No Noise in Self-Reports. \\ ///
	 \end{scriptsize} \end{flushleft} } \end{table} )
eststo clear


erase agg_means.dta
erase cfactual.dta
erase cfactual_results.dta

log close
clear
